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The generic transition in the boson Hubbard model, occurring at an incommensurate chemical 
potential, is studied in the link-current representation using the recently developed directed geomet- 
rical worm algorithm. We find clear evidence for a multi-peak structure in the energy distribution 
for finite lattices, usually indicative of a first order phase transition. However, this multi-peak struc- 
ture is shown to disappear in the thermodynamic limit revealing that the true phase transition is 
second order. These findings cast doubts over the conclusion drawn in a number of previous works 
considering the relevance of disorder at this transition. 
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I. INTRODUCTION 



Quantum phase transitions occurring in bosonic sys- 
tems have experienced a surge of interest lately, due to 
recent experiments showing clear evidence for a quan- 
tum phase transition occurring in an optical lattice^. In 
these experiments, one-, two- and three-dimensional lat- 
tice structures are imposed on the bosonic system by us- 
ing lasers to create standing wave patterns. The lattice 
parameters and the lattice structure of this artificial lat- 
tice can therefore be experimentally tuned, creating an 
almost ideal setting for studying quantum phase transi- 
tions occurring in bosonic systems. Here we shall mainly 
be concerned with the case of a two-dimensional system. 
In the simplest setting the quantum phase transition is 
in this case believed to occur between a Mott insulating 
(MI) phase where the number of bosons per site is fixed 
and the phase therefore incompressible and a superfluid 
(SF) phase where the bosons freely hop throughout the 
lattice. However, several other phases are also possible^. 

In the experiments described above the constituent 
atoms are clearly bosons and quantum phase transitions 
observed in "^He films? should be in the same superfluid- 
insulator (SF-I) universality class. A detailed scaling the- 
ory of the SF-I phase transition arising in two-dimensions 
as described by the bosonic version of the Hubbard model 
was developed in a seminal paper by Fisher et. al^, in 
particular the transition in the presence of disorder was 
investigated and it was pointed out that for the SF-I 
transition in a disordered systems the insulating phase 
would in almost all situations not be a Mott insula- 
tor but instead a compressible Bose glass. Subsequent 
work showed that the same physics should apply also 
to quantum phase transitions occurring in superconduct- 
ing systems dominated by a diverging phase coherence 
length^ making the fcrmionic degrees of freedom irrel- 
evant at the critical point. The reasoning being that, at 
the scale of the diverging phase coherence length the size 
of the individual Cooper pairs would be negligible and 



the underlying physics should therefore be dominated by 
the bosonic degrees of freedom. The scaling theory for 
the SF-I transition^, mainly focusing on phase fluctua- 
tions, should therefore also apply to the two-dimensional 
superconductor-insulator (SC-I) transition. Ensuing ex- 
periments showed support for this scaling picture arising 
in superconducting films^'^'^ and several similar systems 
such as Josephson Junction arrays^*^. This scaling the- 
ory predicts that in two dimensions the quantum phase 
transitions for both the SF-I as well as the SC-I tran- 
sition should take place directly from the superfluid or 
superconducting phase into the insulating phase with no 
intervening metallic phase. However, more recent exper- 
imentsii and theoretical studiesiSiiiiii^ have pointed 
to the possibility of a metallic phase occurring in two- 
dimensional bosonic systems in particular in the presence 
of dissipationiSii^. 

Only at very special points does the SF-I occur at a 
commensurate filling factor. At these points the scaling 
theory of Fisher et al.'^ predicts that the transition is 
in the d + \ dimensional XY class dominated by phase- 
fluctuations. More often the transition occurs at an in- 
commensurate chemical potential which is in a differ- 
ent universality class. This more generic transition is 
expected to be mean-field likc^ and is the focus of the 
present paper. 

Initial quantum Monte Carlo (QMC) simula- 
tionsiSiiLi^ performed directly on the boson Hubbard 
model showed clear evidence for a direct SF-I transition. 
These simulations were constrained to systems with a 
fixed particle number and therefore always incompress- 
ible. Implicitly only the commensurate transition was 
studied. Subsequent studie o^^i^" exploiting a mapping to 
the link-current model removing this constraint focused 
on the phase transition in the presence of strong disorder 
and showed clear evidence for a transition from the 
superfluid to a Bose glass phase, as did studies at fixed 
densities^?. 

A point of controversy has been the phase transition 
occurring at weak disorder. A generalization^ of the Har- 
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ris criterion shows that the transition is stable towards 
disorder \i v > 2/ d. In d = 2 the d + 1 dimensional XY 
does presumably not satisfy this inequality and neither 
does the mean- field v ~ \/2 predicted to occur at the 
generic transition. The scaling theory^ therefore predicts 
that in all cases the Bose-glass phase intervenes between 
the superfluid and insulating phases. This prediction was 
however contradicted by numerical simulations showing 
evidence for a direct SF-I transition at a commensurate 
chemical potential in the presence of weak disorder—. 
Furthermore, Lee et alm^t^^^ reported evidence for a di- 
rect SF-I transition also at incommensurate chemical po- 
tentials. From these studies one would conclude that the 
SF-I transition should be stable towards the introduction 
of disorder, contradicting the scahng theory. 

The recent development of worm algorithmaSSiSi have 
allowed for the study of significantly larger system sizes 
and Prokof'ev and Svistuno\i22iSS have shown evidence 
for the relevance of disorder at the transition occurring 
at commensurate chemical potentials with the universal 
behavior only setting in at very large system sizes . The 
existence of a multicritical line proposed in Ref. I24l25l 
was also ruled out by these studies. However, the con- 
troversy surrounding the incommensurate case still re- 
mains. In fact, to our knowledge, no strong numeri- 
cal evidence exists supporting the evidence of mean-field 
like exponents even in the absence of disorder for the 
generic transition, the results implicit in previous stud- 
ies of the transition in the presence of disorderiSSiSiiS^ 
as well as studies considering longer range interactions 
and the possibility of supersolid phasea^Mi being limited 
to restrictively small lattice sizes. The results including 
longer range interactions were subsequently criticized^^. 
All these studiea^^i^'^i^'*i^^i^^i^^i'^°i'^-^ were performed using 
the link-current representation which we shall also use in 
this study. 

In light of this wealth of experimental and theoretical 
studies it is perhaps surprising that fundamental aspects 
of the SF-I transition as it occurs at an incommensu- 
rate chemical potential are still controversial. This so 
called generic transition is the one most likely to describe 
real experiments and in the present paper we present 
large-scale numerical results using a recently developed 
geometrical worm algorithm, capable of yielding precise 
results for lattice sizes largely surpassing previous stud- 
ies, thereby allowing us to shed new light on this tran- 
sition. In particular we show that, due to the fact that 
this transition is dominated by fluctuations of the par- 
ticle number, simulations on finite lattice will in most 
cases show clearly defined peaks in the energy histograms 
reminiscent of a first-order transition. Only in the ther- 
modynamic limit do the peaks coalesce and a second- 
order transition is recovered. As stressed by Prokof'ev 
and Svistuno v^^i^^ for the commensurate transition, ex- 
treme care should therefore be taken when applying a 
finite-size scaling analysis. For the relative small lattice 
sizes used in the studies by Lee et alMi^ the true influ- 
ence of disorder is therefore likely even further obscured 



by these finite-size energy gaps associated with the first- 
order transition. 

In the remainder of this section we discuss the Boson 
Hubbard model and the particular link-current represen- 
tation that we use for this study as well as the associated 
phase-diagram. Section ^] describes the numerical tech- 
niques employed and Section IIIII focuses on our results 
showing features characteristic of a first-order transition 
for finite lattices, which are found to disappear in the 
thermodynamic limit. We conclude with a discussion in 
Section El 

A. The model 

The simplest model we can write down for interacting 
bosons must at least include an on-site repulsion term 
U as well as a competing hopping term parametrized by 
the hopping strength t. In the absence of the on-site 
repulsion term a bosonic system would always condense 
and would always be superfluid at T = 0. Since we in the 
present paper in particular focus on the generic transition 
occurring at an incommensurate filling we also include a 
chemical potential, fi. We thus arrive at the well known 
boson Hubbard model: 

HhH = J2 (jf'rinr - 1) - ^l^^rj - ^ i^l^r'+c.c) . 

r ^ ^ {r,r') 

(1) 

Here ($r) is the creation (annihilation) operator at 
site r and fir = *&J^r is the number operator. At t ~ Q 
the bosons are completely localized in Mott insulating 
phases while it can be shown"* that a non-zero t eventually 
gives rise to a superfiuid phase with the Mott insulating 
phase persisting in a series of lobes into the superfluid. In 
the MI phase the particle number is fixed and only at the 
tip of the MI lobes does the density not change at the 
SF-I transition. This transition is therefore dominated 
by phase-fluctuations and is in the d+ 1-dimensional XY 
class. The generic transition, occurring at incommensu- 
rate filling factors, is however dominated by fluctuations 
in the particle number and is expected to be character- 
ized by mean-field exponents^. 

We can simplify the Hamiltonian Eq. Q by integrating 

out amplitude fluctuations. We flrst set = l^rle*^"', 
where 9^ is the phase of a quantum rotor. By performing 
the intcgration-'^^'^" i/bH then becomes equivalent to the 
(N=2) model of quantum rotors (cos(6'r), sin(0r)), that 
describes a wide range of phase transitions dominated by 
phase- fluctuations: 




Here, t is the renormalized hopping strength and 7^- = 
Lr is the angular momentum of the quantum rotor. The 
angular momentum can be thought of as describing the 
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deviation of the particle number from its mean, ~ 
rir — riQ. Hence, an equivalent Josephson junction array 
form of this Hamiltonian is: 



r r 

-t C0s(6'r -9r')- 

{r,r'> 



(3) 



If Hqr is written in its Villain- - form we obtain a classi- 
cal modeliSiSSiSliS^ in the same universality class where 
the Hamiltonian is written in terms of integer currents 
defined on the links of a lattice, J = ( J^, J^, J'^) that we 
shall use in this study. One findsi^iSS: 



H 



-E 

K ^ 

(r,r) 



ij2 



r.r) 



(4) 



In this (2+1) dimensional classical Hamiltonian the link- 
current variables describe the total "relativistic" bosonic 
current which has to be conserved on the space-time 
lattice and therefore has to be divergenceless, V • J = 
0. The link-current variables take on integer values 



r^jv,.r = 0,±1,±2,±3. 



Intuitively J: 



(r,T)' 



describe the integer number of bosons hopping in the x 
or y direction from the site r at imaginary time r where 
as Jj"j. denote the number of bosons that remain at the 
site r at imaginary time r. K is the effective tempera- 
ture, varying like t/U in the quantum rotor model. 

In this representation, when K = there is an integer 
number no = (J^) of bosons on each site in order to 
minimize the energy per site. All link variables of the 
classical 3D model vanish in the space directions. The 
"ground state" is composed of no bosons per site when 
the chemical potential is in the interval no — 1/2 < /i < 
no -|- 1/2, and the compressibility k, defined as k 
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zero. This is the incompressible Mott Insulating phase 
with precisely no bosons per site. 

In the other limit K — > oo, the bosons are free to hop 
on the lattice and condensate in the lowest energy mode 
fc = 0. We have off-diagonal long range order and the 
system is in a superfluid phase. The compressibility k, is 
non-zero since the boson number is fluctuating. 



B. Phase Diagram for the pure case 

The phase-diagram of Eq. (Q) has been obtained using 
several equivalent approaches. Sheshadri et alm^ decou- 
ple the hopping term in the following manner: 

$|:$r' = {H)^r' + {^r')H - {H){^r')- (5) 

Using this decoupling and writing ip* — {^l),t/j — (<I>r} 
we see that Eq. can be written as follows: 
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FIG. 1: Mean field phase diagram in the {fi/U, Zt/U) plane 
for the boson Hubbard model 
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(6) 



Here Z is the coordination number. In the number basis 
is non-diagonal but straightforward to diagonalize 
out to fairly large occupation numbers. In this trun- 
cated basis the ground-state energy can be determined 
and the optimal value of ^ determined where the en- 
ergy is minimized. This procedure explicitly yields the 
ground-state wave-function at mean-field level and hence 
directly the density. The superfluid density can be de- 
termined from ps — IV'P- This approach is equivalent 
to preceding work by Fisher et al.^ who showed that the 
Mott insulating lobes with occupation n could be deter- 
mined by considering an action of the typical Landau 
form 5oo(V') = (3N^r\ip\'^ + . . with r given explicitly at 

r = oby 



^ f n+1 n 
r = 2 2 — H ^ — 

Zt \n-fiU -{n-l)+fiU 



(7) 



The transition from the Mott insulator with occupation 
n > to the superfluid is in this mean fleld approach 
given by r = with solution—: 



U 



1x1 



2 2 



± - Vl - 2x(2n+ 1) + 



(8) 



Here, x — Zt/U. For the Mott insulator with n = the 
transition is simply given by ^/U — —Zt/U. Note that 
in Eq. the coefficient in front of — J2i /i/J7-f 1/2. 
Hence the Mott insulating lobes are not centered around 
their corresponding chemical potential, but are off set by 
1/2. The resulting phase diagram is shown in Fig.^and 
is in quite good agreement with detailed strong-coupling 
expansion322i2&. In Fig.^Eq. © are shown along with 
detailed calculations using Eq. ((HJ using the approach of 
Ref.lM 

It is now straightforward to consider the Mott insulat- 
ing lobes present in the quantum rotor model, Eq. |(2Jl, 
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FIG. 2: Schematic phase diagram in the {^/U,K) plane for 
the link-current Hamiltonian (|1J. We have simulated the tran- 
sition between the Mott insulator and the SuperFluid phase 
at /i = 1/4 (arrow in the figure). 



Eq. (PJ. This can be done simply by studying the limit 
n — > oo of Eq. considering that nZt/U remains finite. 
One immediately finds that the lobes are given by: 



U 2 2 



1-4- 



, nZt 
IT 



(9) 



For the link-current model Eq. JJl the only non-trivial 
interaction comes through the global divergence-less con- 
straint and a mean-field treatment is less straight for- 
ward. However, as outlined in Ref. the coupling in 
the link-current model Eq. ^ is related to t/U of the 
quantum rotor model Eq. ^ and the overall shape of 
the phase-diagram must therefore be the same. Since 
the particle number now describes the deviation from 
the mean the resulting phase-diagram therefore becomes 
completely periodic in fi/U and is schematically indi- 
cated in Fig. 121 

In a previous study^^, we studied the quantum phase 
transition at the tip of the lobe fi — Q (black dot in 
figure 10)), where it is known that this model is in the 
universality class of the 3D XY model^i^S. We gave a very 
precise estimate of the critical point Kc = 0.33305(5) 
and found a value — 0.670(3) for the correlation length 
critical exponent, in perfect agreement with the 3D XY 
universality class^. The dynamical critical exponent z 
is equal to unity in this case. 

In this work, we will focus on the pure case where the 
chemical potential is the same for all sites /ir = fJ-- In 
particular, we will present results obtained for a value 
/Lt = 1/4 of the chemical potential indicated by the arrow 
in Fig At this point, the quantum phase transi- 

tion is expected from scaling theory^ to have a dynam- 
ical critical exponent z = 2, and mean field values for 
other exponents (in particular i' = 1/2). As we shall see, 
our simulations confirm this with a good accuracy in the 
thermodynamic limit, but we find finite size effects that 



are strongly reminiscent of a first-order phase transition. 

This has important implications for studies considering 
the relevance of disorder at this generic transition^iSiS^ 
since our results show that for the lattice sizes used in 
these studies the phase-transition appears first-order like. 



II. NUMERICAL METHOD AND 
SIMULATIONS 

A. Numerical method 

We perform Quantum Monte Carlo simulations of the 
model (@J with the recently proposed geometrical worm 
algorithmSLii in its "directed" version^. We refer the 
interested reader to these references for more detailed 
explanations on this numerical technique. We also note 
that our algorithms are closely related to the classical 
worm algorithm by Prokof'ev and Svistuno\i2& (a com- 
parison between the two approaches has been made in 
Ref4i). The most important point to underline is that 
these non local Monte Carlo algorithms permit the study 
of much larger systems with much higher precision than 
what was previously possible using local update schemes, 
thereby getting closer to the thermodynamic limit and 
allowing for much more precise estimates of the critical 
properties of the model. 



B. Simulations 

We are considering a quantum phase transition with 
two different correlation lengths: ^ is the correlation 
length in spatial directions {x,y) and the correlation 
length in the r (imaginary-time) direction. Usually one 
defines ~ thereby defining the dynamical critical 
exponent z, not necessarily equal to unity. To respect 
this anisotropy between space and time directions, it is 
necessary to simulate the the model (0J) on lattices of 
sizes LxLxLt. Periodic boundary conditions are here as- 
sumed and the length of the lattice in the r- direction 
has to be chosen such that Lr = aL^, with a being the 
aspect ratio constant throughout the simulations. This 
follows from the fact that any finite-size scaling function 
will be a function of two arguments: 



f{^/L,^r/Lr)=g{aL,a). 



(10) 



It is therefore necessary to keep the aspect ratio a con- 
stant in order to observe scaling. The value of the dy- 
namical critical exponent is a priori unknown, and one 
usually has to try several values for z to check the validity 
of theoretical predictions. 

Most of the data for different values of the effective 
temperature K were obtained with reweighting tech- 
niques^ by large runs (of the order of 5 x 10^ worms) at 
a single value of K although we in some cases combined 
simulations at several values of K using multi-histogram 
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techniques. We checked every time that our data were in 
the range of vahdity of the reweighting. The error bars 
are obtained with standard jackknife resamphng tech- 
niques^. 

Among several possible thermodynamic variables, we 
have calculated two quantities capable of distinguishing 
the different phases of the systemi2i22. The first quan- 
tity is the stiffness p which characterizes the response of 
the system to a twist in the boundary condition in the 
real space direction. In terms of the link variables, the 
stiffness is calculated a a^^i^°i^^ : 



1 / 2\ 

p = j^M 



(11) 



where — r '^fr t) ^'^^ winding number in the x 

direction. 

The second quantity is the compressibility k, charac- 
terizing the fluctuation of the number of bosons in the 
system. This quantity is written asiSiSS: 



with 



(12) 



being the winding number in 



the time direction. This last quantity can be interpreted 
as the "number of bosons" in the systems. Please note 
that {n^) is in general non-zero when we consider a non- 
zero chemical potential. 

We will use finite-size scaling relation for these two 
quantities to localize critical points and characterize 
them. In two dimensions and near a critical point Xc, 
the scaling theory^ti^iSS predicts the following finite size 
scaling forms for the stiffness and the compressibility : 



L 



2-d-z ~ 



K = L''-'^fi{L^/''S,a), 



(13) 



(14) 



where 6 — \K — Kc\ is the distance to the critical point, 
the correlation length critical exponent and where p and 
k are scaling functions. 

From these scaling forms, we see that right at the crit- 
ical point Kc, in two dimensions d = 2, the quantities 
pL^ and kL'^~^ must be independent of the system size 
if we choose the aspect ratio a to be constant for differ- 
ent lattice sizes. Thus a plot of pL^ or kL^~^ versus K 
for different system sizes should show a crossing of the 
different curves at a single transition temperature Kc- 

Moreover, by simply differentiating the equation H13|) 
with respect to the coupling K, we easily see that at 
the critical point Kc, we have (keeping the aspect ratio 
constant) 



dK 



(15) 



which is used to get the correlation length exponent. We 
find this way of determining v much preferable to the 



traditional data collapso^^i^'^i^^i^^i^^i^^'^^ , which leads to 
significantly more uncertainty in the determination of the 
critical exponents. 

The derivative of p with respect to the effective tem- 
perature K can be obtained by a numerical derivation 
of the curve p{K) or more preferably by calculating the 
thermodynamic derivative during the Monte Carlo sim- 
ulations : 



dp 
dK 



{pE)-{p){E) 
where E is the total energy. 



(16) 



III. RESULTS 

Throughout the rest of this work, we will show re- 
sults obtained with the "directed" geometric worm al- 
gorithmSLii for the link-current model Q at an incom- 
mensurate chemical potential, p = 1/4. We exclusively 
consider the two dimensional case with no disorder. 



A. Determination of critical point 

First, we address the question of the value of the crit- 
ical point Kc and of critical exponents. The scaling 
theory^ predicts z = 2 for the dynamical critical ex- 
ponent. A plot of pL^ or K, should then show a cross- 
ing of the curves for different systems sizes at a single 
point. We have calculated this quantity near the previ- 
ous estimate^SiSi of the transition point Kc — 0.283(3) 
for different lattice sizes for different values of the aspect 
ratio : a = 1/32, 1/16, 1/8, 1/4, 1/2, 1. The small aspect 
ratios allows us to treat larger systems in the real space 
directions, gaining in precision. For all values of a, we 
simulated large systems, up to more than 3 x 10^ lattice 
variables. The maximum size used for different aspect 
ratios is indicated in table ITlL^ A more thorough dis- 
cussion of the infiuence of the aspect ratio follows in sec- 
tion llll 01 Previous work'^°-'^^ using a local algorithm was 
limited to a maximum size L = 10 with an aspect ratio 
a = 1/4 (that is to say lattices of size 10x10x25), whereas 
with the help of the geometrical worm algorithm we have 
been able to simulate lattices of size up to 88x88x242 for 
the smallest aspect ratio used a = 1/32 and in principle 
even larger lattice sizes could be studied. However, as we 
shall see, features present in the results obtained using 
these lattice sizes already tell us that the extrapolation 
to the thermodynamic limit will be difficult. 

In Figure El we show results for the dependence of pL^ 
on K for different values of L and the set of different 
aspect ratios. We see a very good crossing in all cases, 
and the values of Kc (listed in table IIII A|) only show a 
very slight variation with the aspect ratio, and all con- 
verge to give an estimate of Kc — 0.28299(2). As one 
would expect, the universal value of pL^ at the critical 
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a=1/32 a=1/16 




0.28292 0.28295 0.28298 0.28301 0.28304 0.28291 0.28294 0.28297 0.283 0.28303 

FIG. 3: The scaled stiffness pL^ (vertical axis) versus coupling constant K (horizontal axis) for different values of aspect ratios 
a and system sizes L. Lines are guides to the eye. In all cases, the curves for different systems sizes cross at the critical point, 
giving a precise estimate of K^- Insets: at Kc as a function of system size L in log-log scale. Solid lines are power-law 

fits. 



point depends on the aspect ratio (see formula 1 13|l . and 
is listed in table IIIl'XI 

In figure 01 we show our results for the compressibil- 
ity K, as a function of the effective temperature K, for 
all aspect ratios simulated. Also in this case do we ob- 
serve a well defined crossing of curves for different sys- 
tems sizes as one would expect from finite-size scaling 
theory. There is a slight variation of the estimated Kc for 
the different values of the aspect ratios (see table IIII A|l , 
but we can safely estimate Kc from the crossing of k to be 
Kc = 0.28298(2). This estimate of Kc agree within error 
bars with the one previously obtained from the cross- 
ing of pL'^, even if they do not perfectly coincide. It is 
natural to expect such tiny deviations, and only in the 
thermodynamic limit should all estimates (for all differ- 
ent aspect ratios and for all different estimators of Kc) 
give the same value for Kc- The precision at which we are 
able to calculate Kc largely exceeds previous studies^S*^. 

The fact that all our results in Fig Q and Fig. show 
a single well defined crossing for the large system sizes 
used lends strong support support to a value of z = 2, 



as predicted by the scaling theory'*. This value of z was 
implicitly used in the simulations through the choice of 
the lattice sizes and the well defined crossings implies 
that this choice was correct. In the next section (jlll B|l 
of the paper, we will show more convincing numerical 
evidence for this value of the dynamical critical exponent. 

Another important conclusion can be drawn from the 
Fig (O and Fig. (gj: The fact that the estimates of Kc 
from the scaling of pL^ and k show a single well-defined 
Kc allows us to rule out a very hypothetic scenario con- 
sisting of two separate transitions with an intervening 
exotic phase. Our results are clearly only consistent with 
a single well-defined transition. 

We now turn to a discussion of our estimate of v the 
correlation length exponent. In the inset of each of the 
six panels in figure|3| is shown, on a log-log scale, the size 
dependence of calculated at the associated critical 

point Kc- From a power law fit (see equation H15II ). we 
can estimate the correlation length critical exponent v 
for all aspect ratios. The resulting estimates are also 
listed in table IIIl'Xl In the present study we are able to 
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a=1/32 



a=1/16 




0.16 
0.2829 



0.18' 

0.28293 0.28296 0.28299 0.28302 0.28305 0.2829 0.28293 0.28296 0.28299 0.28302 0.28305 



a=1/8 



a=1/4 



0.34 




0.26 



0.22 ' 
0.28292 



0.49 F 



0.41 



0.33 



0.25 
0.28292 



0.28295 



0.28298 

a=1/2 



0.28301 




0.22 

0.28304 0.28292 



0.55 



0.28296 



0.283 



0.28304 



a=1 



0.28296 



0.283 



0.28304 



0.45 



0.35 



0.25 




0.2829 0.28293 0.28296 0.28299 0.28302 0.28305 



FIG. 4: Compressibility k (vertical axis) versus coupling constant K (horizontal axis) for diflterent values of aspect ratio a 
and system sizes L. Lines are guides to the eye. For all aspect ratios, a well-defined crossing of k for different system sizes is 
apparent, giving a value of Kc agreeing with the estimate obtained from the crossing of pL^ (figure EJ. 



simulate significantly larger systems than was previously 
possible in particular for small values of a. For a — 1/32, 
we obtain a value ^ = 0.52(2). This estimate is in very 
good agreement with the mean field value v — 1/2 given 
by the scaling theory of Fisher et a/.'s theoryi. The small 
deviation observed for larger aspect ratios could be due 
to small finite size deviations as we shall discuss in the 
subsequent section. 

To conclude this section, we would like to stress once 
more the importance of simulating larger systems to get 
more precise estimates of the critical properties. A very 
slow approach to the thermodynamic limit is seen for all 
the data in Fig. |j2Jl and Q). For reasons of clarity we 
have not included the data for smaller sizes in these two 
figures. However, the smallest lattice sizes show a clear 
deviation from scaling and the lattice needed to obtain a 
well defined crossing of the curves is in most cases con- 
siderable. In section IlII CI we will show that significant 
systematic finite size effects could easily give rise to an 



incorrect interpretation of data for certain aspect ratios. 



B. Scaling plot for z 

In this section, we will show further numerical evi- 
dence for a value of z = 2 for the dynamical exponent. 
In the context of quantum phase transitions in quan- 
tum spin glasses, Huse and coworkersii (see also Refi^) 
used a numerical method exploiting the anisotropy in the 
imaginary time direction to obtain an estimate for the 
dynamical exponent and the critical point. They pro- 
posed to plot the Binder cumulant versus the aspect ratio 
a = Lt/L^ for different lattice sizes L (i.e. they simu- 
lated systems of size LxLxLr for different Lt). On the 
basis of scaling arguments for the Binder cumulanli^, it 
is then possible to show that all data should collapse onto 
a single curve at the critical point if the correct value of z 
is used to calculate a from the L^- used in the simulation. 
Alternatively, if one for several different values of L can 
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Aspect ratio a 


1/32 


1/16 


1/8 


1/4 


1/2 


1 


Maximum system size Lmax 


88 


64 


56 


50 


40 


32 


Kc (estimated from pL^) 


0.282982(5) 


0.282993(7) 


0.283000(6) 


0.282997(7) 


0.282999(6) 


0.282998(7) 


pL"^ at Kc 


0.397(17) 


0.863(70) 


1.177(60) 


1.177(60) 


0.867(40) 


0.625(25) 


V 


0.52(2) 


0.46(2) 


0.47(2) 


0.47(2) 


0.47(15) 


0.47(15) 


Kc (estimated from k) 


0.282975(8) 


0.282980(12) 


0.28300(12) 


0.282991(8) 


0.282996(8) 


0.292987(8) 


K at Kc 


0.196(6) 


0.227(9) 


0.299(16) 


0.328(10) 


0.374(16) 


0.415(15) 



TABLE I: Estimates of several variables for different aspect ratios used in the simulation. See text for details. 
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FIG. 5: Scaling plot of pL^ versus aspect ratio a = -^j- for 
different system sizes at A" = 0.28299. 



locate a maximal value of the Binder cumulant as a func- 
tion of Lt, then the L™"^ associated with this maximum 
should scale with L as ^ U- . For this to work a 

very precise estimate of the critical point is presumably 
not necessary. 

We adapt this method to the quantum phase transi- 
tion studied here but use the quantity pL^ instead of the 
Binder cumulant. This quantity displays the same prop- 
erties as the Binder cumulant for the application of the 
Huse method but suffers from the drawback that it de- 
pends on an initial assumption of z. Please note that, 
the a priori unknown value of z enters implicitly in both 
axes in this plot. Since we have a good estimate for that 
z = 2, we first calculated pL^ for different L and Lr, 
instead of trying less probable values for the dynamical 
exponent. 

We show the results of our calculations in figure [S] for 
a value of the effective temperature equal to our previ- 
ous estimate of the critical point Kc — 0.28299(2). Please 
note that due to computer time restrictions, we only sim- 
ulated systems of moderate size {L < 40) for aspect ratios 
a < 1/2. Clearly all curves start to collapse into a single 
one, giving strong evidence that our previous estimate for 
z = 2 was correct. This also confirms the validity of our 
previous estimate of Kc- Some deviation from a perfect 



collapse behavior is however present and is attributed to 
small finite size effects since for this calculation we only 
considered systems with size L < 40. In particular, small 
differences in the finite-size estimate of Kc (as was shown 
in table ITlL^ for different aspect ratios are most prob- 
ably at the origin of this small deviation. We further 
discuss finite size effects in the following section. 



C. Effects Depending on the Aspect Ratio 

1. Multiple Peaks 

From the results presented in the previous two sections 
it is clear that the approach to the thermodynamic limit 
is quite slow and pronounced finite-size effects are present 
at small lattice sizes. In the presence of long-range in- 
teractions it is known that the transition in some cases 
can be first-order—. In light of this result we verified 
whether the present transition also showed signatures of 
a first-order transition even though only the on-site re- 
pulsion term U is included. We study this by examining 
in detail the energy distribution close to Kc- In all cases 
we take z = 2. We examined carefully the energy distri- 
bution near the critical point for the six different values 
of the aspect ratio a = 1, 1/2, 1/4, 1/8, 1/16, 1/32 used in 
our simulations. Our results are shown in Fig. where 
we show the probability P{e) of observing an energy per 
site e versus e at the critical point K = Kc- As is evident 
from the results in Fig. multiple peaks are present in 
the energy distribution for the three largest aspect ratios 
while the distribution of the energy for the three smaller 
aspect ratios show a single peak at the critical point. 

Given this observation, it is instructive (see Fig. O to 
extract the contribution to P{e) from the different sectors 
of winding number Ur in the geometrical worm algorithm, 
which can be identified with the particle number (the 
number of bosons) in the system. From the results shown 
in Fig. it is clear that for the three smallest aspect 
ratios, many sectors of the particle number contribute to 
P(e) aX K — Kc, in particular so for the smallest aspect 
ratio a — 1/32. However for the larger aspect ratios, a — 
1/4,1/2,1, only the sectors Ut = 0,1 and eventually 2 
contribute significantly to P{e) at K = Kc- In particular, 
the two main peaks observed in the total histogram P{e) 
can be clearly attributed to the contributions from the 
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sectors with or 1 particle in the system. Moreover, 
these two peaks are clearly separated. 

The generic transition is clearly driven by density fluc- 
tuations as stressed by Fisher et al.^ . If this transition 
is second order, many particle sectors should contribute 
to P{e) at the critical point. This is clearly the case 
for a — 1/32,1/16 and to a certain extent for a = 1/8. 
However, the multi-peak structure in P(e) observed for 
a — 1/4, 1/2, 1 could be interpreted as a signature of a 
first order phase transition. In particular, the energy 
gap between the contribution from the two main parti- 
cle sectors (which is the energy difference between the 
maxima of the two main peaks) should remain non-zero 
in the thermodynamic limit for a first order transition. 
Although it would seem only a remote possibility that 
the order of the phase transition could change with the 
aspect ratio, a more thorough analysis of the signatures 
of a first-order transition observed for a — 1/4, 1/2, 1 is 
clearly needed, as done in the next section. 

2. Scaling of Peaks 

It is not possible to draw any definitive conclusion 
about the order of the phase-transition from the energy 
distribution P(e) calculated for a single L for the three 
aspect ratios a = 1/4, 1/2, 1. It is necessary to verify that 
this is not simply a finite size effect that will disappear 
as the thermodynamic limit is approached. Indeed, since 
we are doing simulations on a finite lattice, the density of 
bosons will not fluctuate continuously but only to vary 
it in steps of There have been many situations 

where one observes evidences for a first order transition 
for small lattice sizes which turns out to be second order 
when simulating larger systems (see for example Ref.-^^). 

We therefore made simulations for several lattices sizes 
for the aspect ratios a = 1/4, 1/2 and 1. We observe for 
all cases a multi-peak structure in the probability distri- 
bution P[e). To treat correctly the size dependence of 
these two peaks, one usually employs the method pro- 
posed by Lee and Kosterlitafi to distinguish between a 
first or a second order phase transition. First, with the 
help of reweighting techniques^, we locate the temper- 
ature K where the two peaks are almost of the same 
heighli^. Then we calculate at this temperature the free 
energy-like quantity — InP(e). 

Our results are shown in Fig. [7| For all the lattice 
sizes L, we still have two pronounced peaks which seem 
to stay of about the same size (or actually slightly in- 
crease for a = 1/4, this is due in this specific case to the 
proximity of the n^. = 2 peak, see footnote*^ and below) 
as L increases. More quantitatively, we plotted in the in- 
set of figure [71 the variation with system size of the "free 
energy difference"— AF — ln( p„,;„|^j ), where P™'*'^(e) is 

the height of both peaks and P™'"(e) the height of the 
minimum between them. We see that this quantity is 
constant with L within error bars. In Lee-Kosterlitz's 
method'*^ , this indicates a second order phase transition. 



To demonstrate this even more clearly, we also note 
that the separation of the peaks (the energy gap) de- 
creases with the system size, as can be seen in Fig. |H1 
In particular, we observe that whereas the first peak 
(corresponding to n,- — 0) stays peaked around a value 
e ^ 0.015 , the peak corresponding to n-r = 1 is clearly 
shifted towards the first peak as we increase system size. 
Looking more carefully at the system size dependence of 
the energy gap Aqi between the peak for the rir = 
and Ur = 1 sectors, we observe that it vanishes as 1/L^ 
(see insets of Fig. We also observe the same scaling 
for the other gaps A02 , A03 corresponding to secondary 
peaks. This behavior clearly corresponds to a second 
order phase transition occurring in the thermodynamic 
limit. The scaling of the separation of the peaks is 
presumably simply due to the fact the particle (boson) 
density varies as between the two peaks on the fi- 

nite lattice. In the thermodynamic limit the density can 
vary continuously. In essence, the spacing between the 
peaks is just a reflection of the finite energy cost asso- 
ciated with adding an additional boson. Moreover, we 
observe that the width of the peaks is decreasing with 
system size. 

We conclude that the observed transition is a second 
order phase transition, but with strong finite size effects 
for large aspect ratios reminiscent of a first-order phase 
transition. These effects are presumably simply due to 
the fact that the density of bosons can not be varied 
continuously on a finite lattice and a result of the finite 
energy required to add an extra particle. Interestingly, 
these finite-size effects do not seem to be refiected too 
much in the estimate of the critical point Kc obtained by 
crossing of pL^ or the compressibility k and could eas- 
ily have been missed in previous studies. However, it is 
immediately clear that if only the sector with particle 
number = is contributing significantly at Kc then 
we are effectively studying the transition occurring at a 
fixed particle number and not the generic, incommensu- 
rate transition. The transition at fixed particle number 
is equivalent to the d + 1-XY commensurate transition 
occuring at the tip of the MI lobes. The non-zero chem- 
ical potential is decoupled and is not taken into account 
correctly. The critical exponents calculated are therefore 
likely to be strongly influenced. This result is of signif- 
icant importance for interpreting what was observed in 
previous simulations where only very modest lattice sizes 
were used"^"'"^^, and the above effects therefore likely to 
have been pronounced. 

These finite size effects could also be of importance for 
previous studies of the generic transition in the presence 
of weak disorderiSSiSiSiS^. The situation is perhaps less 
clear here since a strong enough disorder would smear 
the multi-peak structure observed above and enhance 
particle number ffuctuations. However, weaker disorder 
would only broaden the individual peaks slightly and one 
would effectively be studying the influence of disorder on 
the commensurate transition occurring at /u = (albeit 
at flxed particle number). Recent calculation a^^'^^ done 
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FIG. 6: Probability, P(e), (vertical axis) of energy per site e versus e (horizontal axis) observed during the Monte Carlo 
simulations for different aspect ratios ai K — Kc- Here L — 32 for all aspect ratios. Also shown is the contribution of each 
particle sector (see text) to P(e). For the four larger aspect ratios, contributions of the n,- = 0, Tir = 1 and rir = 2 sectors are 
singled out. For reasons of clarity, we did not denote contributions of higher number of particle sectors. For the two smaller 
aspect ratios, the contributions of all particle sectors are also not denoted for clarity reasons. 



with a varian1iS& of the worm algorithm used in this work 
have shown that the situation in this case is much more 
subtle than previously thought, questioning the vahdity 
of older work considering the influence of weak disorder 
at the commensurate transition22i2i2iiS^. In fact, it was 
observed that the universal behavior only sets in at very 
large lattice sizes. Hence, if only weak disorder is consid- 
ered for small lattice sizes one is therefore likely to ob- 
serve incorrect exponents. This would explain why Park 
et alm^ observe z/ = 0.5 ± 0.1 at the generic transition 
in the presence of weak disorder violating the inequality 
V > 2/d and would question the validity of the phase- 
diagram presented by Lee et alM^. 



D. Correlation Functions 

Even though in the previous two sections we have 
seen that the pronounced finite-size effects present at the 
larger lattice sizes have relatively little influence on k 
and pL^ they are quite visible in the correlation func- 
tions. The correlation functions can be easily calculated 
using the geometrical worm algorithm as described in 
reference |4l[ From scaling theor}4i2£ we expect them to 



follow the following form at K — K^. 

a(r) -T-(''+^-2+'')/^ (17) 

In Fig. 1^1 we show representative results for Ct(t) 
calculated dX K = Kc for the generic transition at 
/i — 1/4. There is a clear power-law dependence and 
it is relatively easy to extract the associated exponent 
Ct{t) = 0.47742r~°-^^*. If we take our previous esti- 
mate of z = 2 for granted and remember that d — 2 
this would imply that 77 = consistent with mean-field 
behavior. 

The behavior of the correlation functions in the spa- 
tial direction is much more complicated. In Fig. (|10ll we 
show representative results for Cx calculated at K = Kc 
for a system of size L = 32 for a range of different 
Lr = 16, 32, 64, 128, 256, 512. As is clearly evident from 
this figure these correlation functions do not follow a 
power-law behavior. Given the results presented in the 
previous sections this is not surprising and illustrates 
the difficulties associated with these finite-size effects. 
The deviation from power-law behavior increases dra- 
matically with the aspect ratio. Hence, only by studying 
significantly larger system sizes at much smaller aspect 
ratios would one presumably be able to recover the ex- 
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FIG. 7: The free onergy-likc quantity — InP(e) (sec text) for the three aspect ratios a = 1/4, 1/2, 1 versus energy density e. 
The curves for different system sizes L have been arbitrarily shifted vertically for clarity reasons. Insets : The free energy 
difference AF (see text) versus system size L. 




FIG. 8: Scaling of the peaks in P(e) for different system sizes for the aspect ratios a = 1/4, 1/2, 1. Insets : Energy gap Aoi 
separating the two first peaks as a function of IjL^ . Dashed line is a a fit to Aoi = Ajl? . 
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FIG. 9: The r-dependent correlation function calculated on 
a lattice with L — 64 and Lr — 128 at Kc- The dashed 
line (hidden by the data points) indicates a fit to the form 




X 



FIG. 10: The x-dependent correlation function at A" = Kc for 
L — 32 calculated for a range of different Lr — 16,...,512. 



pected power-law behavior. We have so far been unable 
to do so. 



IV. CONCLUSION 

In conclusion, we have in the present paper presented 
large-scale numerical results for the SF-I transition oc- 
curring in the boson Hubbard model at an incommen- 
surate chemical potential /U = 1/4 obtained using the 
"directed" version of the recently developed geometrical 
worm algorithmSLii. By carefully analyzing the proba- 
bility distribution P(e) of the energy density at K = Kc 
for different values of the aspect ratio, we showed that 
strong finite-size effects, reminiscent of a first order tran- 
sition are present for the larger aspect ratios. These ef- 
fects disappear in the thermodynamic limit and the tran- 
sition is indeed second order with mean field exponents 
as predicted by scaling theory**. However, if only small 
lattice sizes are used then these effects are pronounced 
and would imply that the chemical potential is not taken 
into account correctly, resulting in the associated critical 
exponents being calculated incorrectly. Finally, we note 
that amplitude fluctuations, absent from the Hamilto- 
nian Eq. Q used here, could be important at the generic 
transition as recent studies indicate^. 
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